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ABSTRACT 

Computations of viscous-inviscid interacting internal flowfields are presented for steady and un- 
steady quasi-one-dimens ional (Q1D) test cases. The unsteady Q1D Euler equations are coupled 
with integral boundary-layer equations for unsteady, two-dimensional (planar or axisymmetric), 
turbulent flow over impermeable, adiabatic walls. The coupling methodology differs from that used 
in most techniques reported previously in that the above mentioned equation sets are written as a 
complete system and solved simultaneously; that is, the coupling is carried out directly through the 
equations as opposed to coupling the solutions of the different equation sets. Solutions to the 
coupled system of equations are obtained using both explicit and implicit numerical schemes for 
steady subsonic, steady transonic, and both steady and unsteady supersonic internal flowfields. 
Computed solutions are compared with measurements as well as Navier-Stokes and inverse bound- 
ary-layer methods. An analysis of the eigenvalues of the coefficient matrix associated with the qua- 
si-linear form of the coupled system of equations indicates the presence of complex eigenvalues for 
certain flow conditions. It is concluded that although reasonable solutions can be obtained numeri- 
cally, these complex eigenvalues contribute to the overall difficulty in obtaining numerical solutions 
to the coupled system of equations. 


I. INTRODUCTION 

The study and analysis of internal flows has received significant attention over the past several de- 
cades because the operation of many physical devices, particularly regarding aerospace-related 
hardware, depend upon proper designs to achieve near-optimum operating characteristics. Exam- 
ples of such devices include any configuration where the flow is confined and an exchange between 
pressure and kinetic energy is desired (engine inlets, wind tunnel diffusers, rocket nozzles, etc.). 
There devices can be geometrically complex as well as very viscous-flow dominated. Moreover, 
certain configurations and conditions can result in unsteady flow (e.g., inlet buzz). 

In the past, the design of these devices has, for the large part, depended upon empirically based meth- 
odologies. More recendy, computational techniques have played an increasingly important role in 
the design process as hardware becomes less conservative and is required to operate near the edge 
of the design envelop. As evidenced above, perhaps the most important physical flowfield charac- 
teristics which need to be considered when attempting to computationally address internal flows are 
effects associated with unsteadiness, viscosity, and multi-dimensions. Of course, the relative con- 
tributions of these effects are dependent upon the geometry as well as which physical flowfield pa- 
rameters are required to provide the “answers” for a given problem. For example, if the performance 
(e.g., static pressure rise) of a subsonic axisymmetric diffuser is desired, it is very important that 
viscous effects be well represented because diffuser performance is very sensitive, for example, to 
the incoming blockage caused by the presence of the boundary layer. For this case, it can be argued 
that unsteadiness and multidimensional effects play a secondary role. However, for cases where 
boundary-layer separation is possible, significant unsteadiness may be present. For these cases the 
capability to capture this unsteadiness within computation is important in order to gain engineering 
insight into the physics. On the other hand, it is easy to identify cases where all of the above effects 
play an important part in shaping the overall flowfield structure (e.g., a moving shock within an S- 

shaped, asymmetric duct). 

A thorough computational investigation of flowfields of this type requires solution of the full Re- 
ynolds-averaged, multidimensional, time-dependent, Navier-Stokes equations. Of course, solu- 


tion of these equations produces essentially all pertinent flowfield parameters. Therefore, assuming 
that these solutions are of acceptable accuracy, it possible to perform parametric studies of a pro- 
posed geometry/flowfield combination which could be used to significantly reduce the risk 
associated with new hardware design. Unfortunately, obtaining numerical solutions to these equa- 
tions for complex geometries and unsteady flowfields is expensive and time-consuming, even using 
today’s largest and fastest supercomputers. Therefore, it is important to investigate alternative 
means of performing compute-based parametric studies of proposed new hardware designs. How- 
ever, it is equally important that these alternative techniques be capable of capturing as much of the 
critical physics as possible to avoid “throwing the baby out with the bath water.” Consequently, iden- 
tifying the physical aspects which tend to dominate the behavior of the flowfield associated with a 
particular geometry is vital to the success of the alternative computational procedure. 

It is obvious that some compromises must be made to reduce these computational requirements 
while simultaneously retaining the desired physics. Deciding upon which compromise requires an- 
swering the following question: “What are desired physics?” or stated another way, “Which physical 
flowfield characteristics are we willing to approximate in order to reduce the overall computational 
resource demands?” Unfortunately, the answer to either question is very problem dependent. For 
example, elimination of viscosity effects from the Navier-Stokes equations results in the Euler 
equations. Obviously, this reduced^quation set by itself can never be used to simulate the flow of 
a viscous fluid, but can, however, be used to generate “reasonable” solutions for unsteady flow about 
extremely complex, three-dimensional geometries, as demonstrated by Whitfield, et al 1 who com- 
puted the unsteady flow about three-<iimensional transonic propfans using the Euler equations. The 
precise meaning of “reasonable” relates to the above question(s). That is, an assumption was made 
in Ref. 1 , a priori, that viscous effects could be neglected for the configurations and flow conditions 
to be investigated. Comparisons between measured and computed performance parameters 1 indi- 
cated that this assumption was indeed “reasonably” valid. Therefore, it could be argued that the com- 
promise made to exclude viscous effects from the analysis did not contaminate the computed solu- 
tions to the point of being unusable. However, it should be pointed out that this conclusion is based 
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on the original stipulation that (as an example) effects of viscosity and the ensuing ramifications of 
its presence were of lower priority in the simulation. 

Of the three aforementioned physical characteristics under consideration, the one likely to have the 
most significant impact upon computational resource requirements concerns that of multidimen- 
sions. This can be argued from the standpoint that the number of floating point operations required 
for a given simulation is roughly proportional to (»*°")(2 + ndim) 2 , where n is the number of grid 
points and ndim is the number of spatial dimensions. Of course, this proportionality is greatly depen- 
dent upon the numerical scheme used to solve the equations, but at least gives an indication of how 
quickly the cost of performing multidimensional simulations escalates. Similar to the arguments 
given above for three-dimensional viscous and inviscid flows, the validity of compromise (i.e., re- 
duction) in the number of independent spatial variables is problem dependent and is difficult to 
judge, a priori, whether the resulting simulation adequately represents reality. As stated by Hirsch 2 , 
“In all cases, however, the final word with regard to the validity of a given model is the comparison 
with experimental data or with computations at a higher level of approximation.” Therefore, it is 
the reduction in effects associated with multi-dimensions, while retaining effects of unsteadiness 
and viscosity, and solution of the resulting equations (for internal flows) which forms the basis and 
underlying motivation of the present effort. In particular, the development of an engineering tool 
through which preliminary estimates of unsteady internal flow processes can be generated using 
available workstation-based hardware is sought. 

One approach to achieve this is to seek solutions to the unsteady, two-dimensional Navier-Stokes 
equations, or the unsteady two-dimensional Euler equations coupled with the steady (or unsteady), 
two-dimensional boundary-layer equations. While these are valid approaches, even the two-di- 
mensional equations can result in nontrivial computational time requirements, particularly for un- 
steady flow. However, use of the coupling approach (e.g., Euler coupled with boundary layer) has 
significant resource-saving advantages over that associated with solving the full Navier-Stokes 
equations because of relaxed grid requirements in viscous regions 3 . Hence, the coupling approach 
is adopted here, where equation sets valid for a particular region of the flowfield are used. Specifi- 
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cally, the Euler equations written for unsteady, quasi-one-dimensional (Q1D) flow are coupled with 
integral boundary-layer equations for unsteady, two-dimensional turbulent flow over adiabatic 
walls. The assumption is made that solutions to the coupled equations will yield results of engineer- 
ing accuracy. It must be emphasised 4 that the validity of using the simplified equations is very prob- 
lem dependent and, similar to other analytical or computational techniques, requires experience and 
engineering judgement with regard to whether the approach and/or computed solutions represent 
reality. No attempts are made at quantifying specific classes of problems for which the approach 
presented herein can be used. Attempts are made to quantify the validity of these assumptions (or 
lack thereof) through comparisons with available experimental and computational sources. 

An additional assumption fundamental to the coupling approach applied to internal flows is that the 
flowfield within the channel contains an inviscid “core” (i.e., not fully developed) of fluid which 
is allowed to interact with the viscous region near the wall. A schematic of this type configuration 
is shown in Fig. 1. The displacement of mass caused by the presence of this viscous region has a 

thickness of 6*, defined by 


* 

QeUed = 


c 00 

Jo 


(( QeUe - QU) <ty 


which is exact for planar flow, but is only approximate for the axisymmetric case. However, the 
above expression approaches the true mass defect length for axisymmetric flow when the local 
boundary layer is thin compared to the local body radius 5 ’ 6 . Therefore, the analysis presented herein 
is valid only for those cases where the boundary layer is small relative to the local body radius. 


Results ensuing from this analysis is reported here in two parts. Part 1 is essentially a continuation 
of efforts reported in Refs. 4 and 7 where an explicit numerical scheme was used to solve the system 
of equations formed by writing the viscous and inviscid equations as one complete system. This 
coupling methodology differs from those reported previously (e.g., 3 ’ 8 ’ 9 ) where the coupling was 
performed between the solutions to the equation sets rather than the equations themselves. As dis- 
cussed in Refs. 4 and 7, this approach is motivated by the observation that coupling the solutions 
results in a scheme which can have convergence difficulties and is often not robust, particularly for 
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“strong” interaction cases. In addition, previous coupling schemes which use the steady, direct form 
of the boundary-layer equations to solve for the viscous region for cases where boundary-layer sep- 
aration occurs fail because this form of the boundary-layer equations are singular at or near separa- 
tion 3 . (It should be noted though that the singularity can be avoided by using the so-called inverse 
form of the equations 3 . However, because the formulation of the unsteady inverse form is not 
unique, coupling of the viscous and inviscid equation sets is less than straight-forward 10 ). As shown 
by Moses, et al n , however, a simultaneous solution procedure (using the steady form of the laminar 
boundary-layer equations and Laplace equation for the stream function) apparently removes the 
separation singularity which makes the computation of separated flows possible using the direct 
form of the boundary-layer equations. In an analogous manner, the present approach simultaneous- 
ly solves the unsteady forms of the Q1D Euler equations and the unsteady integral boundary-layer 
equations for turbulent flow for steady subsonic (both separated and attached) flows, as well as un- 
steady supersonic (attached) flow cases. However, as discussed in subsequent sections, several of 
the disadvantages which the present direct coupling approach sought to overcome have been re- 
placed with other, perhaps more disheartening ones with regard to seeking numerical solutions of 

the complete system of equations. 
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II. PART 1 SUBSONIC AND SUPERSONIC FLOWS 

a. Formulation of Equations 

Much of the following material is given in Refs. 4 and 7 but is repeated here to provide the necessary 
background for the analysis presented in Part 2, and also for convenience to the reader not familiar 

with previous efforts 4,7 

The equations which form the basis of the present analysis are the unsteady momentum and mean- 
flow kinetic energy integral boundary-layer equations for turbulent flow (for the viscous region, 
Eqs. (1) and (2)) and the unsteady Q ID Euler equations for no work and adiabatic flow (for the mvis- 
cid region, Eqs. (3), (4), and (5)), and are written as a complete system as 4,7 : 


■^(QeUed*) 

- Uej-^eSg) + 

+ = 0 

(1) 


Ue(0 + d* ~ 0p)] 

+ 

2 QeUeOg ~ 


+ 2 Q e Ue(d 

* 

dl)^ - leeuf-E = 0 

(2) 



+ ^(QeUeA) = 0 

(3) 


jjiQeUeA) + 

d 

dX 

(Geul + Pe)*] ~ Pe = 0 

(4) 


+ i 

^[(u e (E e + p e )A] + Pe ^ ~ 0 

(5) 


where A = area, q = density, u = velocity, p = pressure, E = total energy, t = time, and x = axial dis- 
tance. Subscript e has been added to the gas dynamic variables to denote “edge” values, taken to 
be those associated with the inviscidcore. It should be noted that the area A, by definition, represents 
that part of the flow region which contains inviscid fluid. These equations have been non-dimen- 

sionalized using the dimensional parameters /^(length), <? 0 ,« (density), and a 0 co (velocity) where 
a is the speed of sound and subscript 0, « indicates upstream stagnation conditions, and A denotes 
a dimensional quantity. In Eqs. (1) and (2), k = 0 for planar flow or 1 for axisymmetric flow, where 
the following integral length definitions have been used. 
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(6a) 


(5* 



0 = 



0* 



(6b) 

(6c) 


K = \ (i - ft) 

Jo 



The skin friction and dissipation integral are given by 

2t w 

Cf ~ QeUj 

D - f ’*$(*)* 

Jo 

where t and t w are the local and wall shear stress, respectively. 


(6d) 

(6e) 


(7) 

( 8 ) 


To place the system of equations in a form amenable to numerical solution using explicit schemes, 
the temporal derivatives of Eqs. (1) through (5) are isolated, i.e.. 


jt(Q eUe d*) - u e j- t (Q^o) ~ b \ 

(9) 

4:[(eeul(0 + 6* - 0p)] - 2 QeUeiOg ~ = b 2 

(10) 

II 

3 

(ID 

jjiQeUeA) = b 4 

(12) 


(13) 


where bi through b 5 are defined by referring to Eqs. (1) through (5), respectively. 

Reducing the equations further requires choosing a dependent variable vector. In Refs. 4 and 7, the 
dependent variables used were 
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q = (, Qe u e Me 9 H? < 14) 

Once this choice is made, the temporal derivatives are expanded and the ensuing terms algebraically 

manipulated to form a system of five simultaneous partial differential equations, which can be writ- 
ten as: 


where L is a 5 x 5 matrix, and b is the right-hand-side vector containing spatial derivatives. Here 
M e is the edge Mach number, 9 is the momentum thickness (Eq. 6b), and H is a shape factor defined 
using integral lengths formed with kinematic properties, i.e.. 



where 


(16) 



(17a) 

(17b) 


The number of unknown parameters in Eq. (15) is ten, including those contained in the right-hand- 


side vector b. Closure of the system of equations requires that all variables be expressed in terms 
of the dependent variables. This was partially accomplished in Ref. 7 (and here as well) through 
the use of several auxiliary relations involving boundary-layer integral length shape factors and the 
perfect gas equation of state. Complete closure was accomplished in Ref. 7 by using the relation 


^ = 1 + 7 -^-m\ ( lg ) 

T e 1 

which is justified so long as the gas is thermally perfect and that the gas is in a state of equilibrium 
at each cross-section. To illustrate how the above relation was used in Ref. 7, consider the following 
expression for the total (stagnation) energy, written in dimensional form as 


* A A 1 A *2 

E e = Q e e e + 

where e e is the internal energy given by (for a perfect gas) 


( 19 ) 
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( 20 ) 


e, = c v T e 


and c v is the constant-volume specific heat. It is obvious that the temperature must be expressed 
in terms of the dependent variables. This was performed in Ref. 7 as follows. In non-dimensional 
form, the internal energy becomes 

A T 'T T 

T e T 0,e _ 1 T e l Q,e 


e e 

*2 

a 0 ,« 


CyT e 

a 2 

do, oe 


[yR 


y(y " 


( 21 ) 


00 


where T n is the local stagnation temperature and R is the perfect gas constant. Making the defini- 
tions 


results in 


ip- - i + y ; 1 ^ 

T e 1 

(22a) 

. = [y(y - i)/e] — 1 

(22b) 

A 

^ p « rri 

72 e ‘ r o.« 

(23) 

ao,oo 


sion for total energy can be written 


= Qe(&eTQ e + 2 W *) 

(24) 


where 


T 0,e ~ 


[ 0 t e 


’ 0, GO 


(25) 


Differentiating the above expression for total energy results in derivatives of the quantity To, e which 
must be assumed constant or otherwise specified. Therefore, the assumption was made in Ref. 7 
that the ratio of local to reference stagnation temperature was equal to one (constant), thus eliminat- 
ing these derivatives. In addition, the above ratio also appears in several terms within the elements 
of the matrix L. Although this is a good approximation for steady flows involving no heat transfer, 
the validity of the assumption becomes questionable for unsteady flows, even with no heat transfer. 
This can be seen by rewriting the energy equation (Eq. (5)) in terms of stagnation enthalpy, 
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( 26 ) 


|(<M*o> + fx {e ‘ u < Ah o) = 

where ho is the stagnation enthalpy defined by 

h 0 = he + jul ( 27 ) 

Use of the continuity equation results in 

dh 0 dh 0 i d£e (28) 

dt + Ue dX Qe dt 

We can rewrite the left-hand-side of the above as a material derivative to finally arrive at 


Dhp _ 1 dPe (29) 

Dt Qe dt 

which is valid in the absence of work and heat flux. By using the energy equation in this form, it 
is straight-forward to see that changes in static pressure due to unsteadiness results in corresponding 
changes in the stagnation enthalpy, and thus the stagnation temperature. 

This situation can be avoided in at least two ways. One method is to express the internal energy 
instead as 


_ Cv Te _ 
<2o,» y&t o >00 


y(y - i) 


where 


From the definition of the sonic velocity 


a 2 e 


T e = 


0, 00 


a 2 

a 2 

# 0 , 0 ° 


— e — 


= Te 


0 ,® 


and using the Mach number, it follows that 


e e 

*2 

a o,« 


u. 


y(y - ml 


(30) 


(31) 


(32) 


(33) 


which is the desired result. 
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Another way to avoid making the constant stagnation temperature assumption is to replace Mach 
number with static pressure as a dependent variable. To illustrate this, consider the perfect gas equa- 

tion of state given by 


Pe 



or, in non-dimensional form, 


Pe 


eJe 

y 


Therefore, 



It follows that the Mach number can be computed as 


(34a) 


(34b) 


(34c) 


which again permits closure of the system. 

Results ensuing from numerical schemes based upon all three formulations are presented in this re- 
port. As discussed in subsequent sections, both explicit and implicit numerical methods have been 
implemented. Whereas an explicit scheme has been utilized in all formulations, the implicit method 
has been applied to only the formulation where pressure is used as a dependent variable. 

Because of the possible ramifications regarding solutions computed using the original formulation 
which assumes constant stagnation temperature 7 , a brief diversion will be taken at this point to inves- 
tigate differences between computed steady and unsteady solutions resulting from the different for- 
mulations. As mentioned above, numerical computations involving all formulations have been per- 
formed using an explicit scheme (implementation of this scheme is discussed in Section H.c). To 
investigate differences in computed solutions ensuing from these formulations, results for a constant 
area axisymmetric duct ( 10 radii in length) with a fixed entrance Mach number of 2 and a reference 
Reynolds number of 5 million are presented. Converged (steady-state) solutions from these com- 
putations are shown in Figs. 2 and 3. Figure 2a presents time histories of exit Mach number and it 
can be seen that identical results are obtained at steady state for all formulations. Also, convergence 
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of both formulations which do not assume constant T 0 e are seen to be essentially identical. Although 
not shown here, all other pertinent flowfield parameters converge in a similar manner. In addition. 
Fig. 2b illustrates convergence as measured by the root-mean-square residual of the velocity, de- 
fined as 



where N m « is the total number of grid points. It is of interest to note that although convergence to 
machine epsilon (using double-precision floating-point operations) is achieved in approximately 
40 (non-dimensional) time units (corresponding to approximately 400 - 700 time steps at a CFL 
of 0.9), referring back to Fig. 2a indicates that the solution at the exit has stopped changing appreci- 
ably in less than half that time. 


As expected, the different formulations do not give the same results for unsteady flows, although 
differences are observed not to be large, at least for the cases examined thus far. Computed solutions 
from all formulations for an unsteady flow are illustrated in Fig. 3. This is a contrived test case pres- 
ented in Ref. 7 for the duct mentioned above, where now the entrance Mach number is sinusoidally 
varied at a non-dimensional frequency of 0.1 to yield an entrance Mach number with a mean value 
of 2.0 and an oscillatory magnitude of 0.2. Fig. 3 a compares computed exit Mach number solutions 
based upon all formulations where the entrance (input) Mach number is also shown for comparison. 
It can be seen that solutions from formulations not involving the assumption of T Q = constant are 
identical. Although these solutions are very similar to those from the original formulation, mini- 
mum and maximum values of computed Mach number are seen to vary. Also, there is a very slight 
phase shift regarding the (min,max) values of Mach number. However, as shown in Fig. 3b, differ- 
ences in computed boundary-layer integral lengths are much less. The origin of these relatively 
small differences can be seen in Fig. 3c which illustrates the time variation of stagnation temperature 
at the duct entrance and exit, the former of which is constant, by definition. It can be seen that at 
this axial location, the maximum variation of stagnation temperature is approximately three to four 
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percent Therefore, an assumption of T 0 , e = constant (used in the original formulation) is reasonable, 
at least for this degree of unsteadiness. 

Although differences exist between computed unsteady solutions that are based upon different for- 
mulations of the system of equations, these differences are remarkably small, at least for the test case 
shown. However, because the formulation which assumes T 0 , e = constant is inconsistent with regard 
to the simulation of unsteady flows, either of the new formulations are the preferred methodologies 
in this regime. This is particularly true if the present coupling methodology is to be applied to cases 
involving both unsteady flow and added heat flux. 

As discussed previously, the system of equations is written in terms of the coefficient matrix L, where 
the elements of L vary according to the particular formulation. Elements of the matrix L (as well 
as those of another matrix N, to be discussed next) are given in the Appendix for all formulations 

discussed herein. 

b. Eigenvalue Structure of the System of Equations 

In Refs. 4 and 7, the approach to solve the system of equations (14) was to use semi-discretization 

which results in a system of ordinary differential equations at each mesh point. The equations were 
then solved with a two-stage Runge-Kutta scheme using first-order backward spatial differencing 
throughout the computational domain. The exclusive use of upwinding was possible in Ref. 7 be- 
cause for the cases considered, all eigenvalues of the coefficient matrix L~ x N were found to be 
positive (as well as real, of course). The matrix L~ l N results from writing the system of equations 
(14) in quasi-linear form 

d S + L -x N *l = L-'d 01) 

dt dx 

where the matrix N is derived in similar fashion as is the matrix L. Also, the different formulations 
discussed previously result in differences in various elements of the matrix N (similar to the matrix 
L). Elements of the N matrix for all formulations are given in the Appendix. 

It is of interest to examine the behavior of the eigenvalues over the expected range of the various 
parameters upon which elements of the matrix L ~ X N depend. Using isentropic relations between 
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local static and stagnation conditions, it is straight-forward to show that these elements (and thus 
the eigenvalues) can be expressed in terms of M e , H, and 9. Eigenvalue distributions as a function 
of Mach number are shown (as lines) in Fig. 4 with both H and 6 used as parameters. The eigenva- 
lues shown were computed using the elements of the L and N matrices resulting from the formulation 
where M e was used as a dependent variable but T 0 e was not assumed constant. (It is of interest to 
note that eigenvalues computed using the formulation where p e is a dependent variable are virtually 
identical to those shown). Because of the algebraic complexity of the matrix L ~ l N, the eigne values 
were computed numerically using an iterative technique 12 . Also plotted in these figures (as sym- 
bols) are the eigenvalues associated with the Q1D Euler equations written as an isolated system (Ue , 
Ue +Oe, and u* - a*). It should be noted that Reynolds numbers were evaluated assuming reference 
(stagnation) temperature and pressure to be 520 °R and 14.7 psia, respectively. This results in mo- 
mentum thickness Reynolds numbers ranging from approximately 400 to 480,000 for 0.00 1 < 9<0. 1 . 

Figures 4a-4e (0=0.001) illustrate eigenvalue behavior for 1.2 <> H <6.0. It can be seen that for 
shape factors less than approximately 2.0, all eigenvalues remain positive for Mach numbers greater 
than one thus confirming observations in Ref. 7. However, this is not the case for higher values of 
H as shown in Figs. 4d and 4e which indicate at least one eigenvalue becomes negative for Me > 
1. Also, it is interesting to note that three eigenvalues of the complete system closely approximate 
those of the inviscid equations for all values of 77. 

Perhaps the most interesting (or disturbing) aspect of eigenvalue behavior can be seen in Figs. 4d 
and4e which indicate the appearance of complex conjugate pairs at high shape factors and superson- 
ic Mach numbers, where the range of Mach numbers within which this occurs decreases with in- 
creasing H. While only the real part is plotted, the imaginary part is observed to be at least one order 
of magnitude smaller than the real part. It should be noted that the appearance of complex eigenva- 
lues seem to occur for shape factors high enough to be indicative of boundary-layer separation. Dis- 
cussion regarding the ramifications of the appearance of complex eigenvalues is given at the end of 

this section. 
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One might hope that complex eigenvalues would occur only within a relatively small range of values 
associated with the various parameters. Unfortunately, this is not the case which is illustrated in Figs. 
4f through 4o. InFigs.4f-4j (1.2 < H <6.0, 0=0.01), complex eigenvalues again appear, but only 
for shape factors high enough to cause boundary-layer separation (which generally occurs for 
2.8 <1 H < 3.0, depending upon the Reynolds number). The range of Mach numbers over which 
this occurs is rather extensive at higher values of H . It should also be noted that significant deviation 
from the inviscid eigenvalues has occurred at the higher values of momentum thickness, particularly 
for higher shape factors. In addition, negative eigenvalues occur over the entire range of Mach num- 
ber, again for higher values of H. This behavior is even more pronounced for very high values of 
momentum thickness as shown in Figs. 4k-4o (1.2 < H < 6.0, 0=0.1). However, one could argue 
that for 0=0.1 we have violated one of the fundamental assumptions regarding use of the coupling 
approach; i.e., recalling that 0 has been non-dimensionalized by a reference length (usually the inlet 
radius or half-height), a value of 0=0.1 indicates that the local momentum thickness is 10% of the 
local radius. Assuming a shape factor of 1.5 and that the local radius variation is small compared 
to that of the inlet (which is consistent with the Q1D assumption) implies that the displacement thick- 
ness occupies 15% of the channel radius. Assuming further that the displacement thickness is 
approximately 1/6 of the total viscous region implies that the flow is essentially fully developed 
which, of course, violates our original stipulation that this not be the case. Therefore, Figs. 4k-4o 
should be interpreted as an illustration of how the eigenvalues behave toward the upper end of valid 
parameter space. On the other hand, values of 0 in the range of 0.1 do not necessarily imply that 
the channel is fully developed. For separated flows, integral lengths have the tendency to grow rap- 
idly because of the large, retarded flow region near the wall. However, the overall viscous region 
can remain small enough such that an inviscid core exits. Example of this are illustrated in subse- 
quent sections. 

The appearance of complex eigenvalues indicates that the system of equations in their present form 
cannot allow solutions as a well-posed initial/boundary-value problem by integration over time. 
This conclusion is based upon the work of Briley et al 13 who used the criterion set forth by Garabe- 
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dian 14 “that it is natural to require that every root of the characteristic equation be real, as this ex- 
cludes solutions thatmay grow exponentially with the time-like variable” 13 . Therefore, exponential 
growth in solutions of a system of equations (which supposedly represent certain physics) can be 
attributed to numerical instability of the unsteady solution algorithm and/or to a mathematical set 
of equations which is ill-posed for solutions as an initial value problem in time 15 . Differentiating 
between these two areas of concern for the present effort requires that the problem be separated into 
its individual pieces, namely, physics, mathematics, and numerics. That is, the objective is to obtain 
a valid mathematical representation (equations) of the physics and to numerically solve these equa- 
tions in a stable manner to a specified order of accuracy. With regard to the physics of the inviscid 
flow, it has been well established that the Euler equations represent a very good approximation to 
the motion of a fluid in regions where effects of viscosity are neglible. However, with regard to the 
physics of the boundary-layer flow, Whitfield 16 encountered complex eigenvalues in seeking solu- 
tions to the unsteady integral boundary-layer equations where time was used as an iterative parame- 
ter to reach steady state. In addition, similar eigenvalue behavior was encountered in Refs. 10 and 
17 in dealing with the unsteady, three-dimensional integral boundary-layer equations which, how- 
ever, did not preclude obtaining reasonable numerical solutions 10,17 . Along these same lines, the 
integral boundary-layer equations of the type used herein have been shown to yield good engineer- 
ing approximations to viscous flows in regions where the usual boundary-layer assumptions are val- 
id for both steady and unsteady regimes 6,18 . 

Based upon the above discussion, it is reasonable to conclude that it is the approximate governing 
equations which are the origin of the observed anomalies. While the system of equations used m 
the present effort does not generally exclude solutions exhibiting exponential growth due to ill— 
posedness, it is shown in subsequent sections that numerical solutions of engineering accuracy can 
be obtained for those cases where either: (a) no complex eigenvalues are encountered, or (b) if com- 
plex eigenvalues do appear, unbounded growth can occur but can be very slow, thus allowing reason- 
able solutions to be computed. 
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Because the method presented herein uses many of the same shape factor correlations and auxiliary 
relations employed in Refs. 10 and 16, it is believed that the appearance of complex eigenvalues can 
be attributed to the approximations introduced by these empirical and analytical relations; i.e., these 
empirical relations and approximations are insufficient to define a well-posed set of approximate 
governing equations. Of course, this situation calls for an analysis similar to that performed by 
Briley 13 to attempt to locate the specific relations which cause the observed eigenvalue behavior. 
Unfortunately, time and resource limitations preclude pursuing such an analysis. 

c. Numerical Method 

Based upon the preceding discussion, a numerical scheme utilizing spatial difference operators other 
than purely one-sided is required for the general case unless, of course, a completely upwind method 
which uses spatial differences whose type depends upon local flowfield characteristics (eigenvalues) 
is used (note this approach is not even applicable for situations resulting in complex eigenvalues). 
Because of the algebraic complexity of the governing system of equations and the above concerns 
regarding well-posedness, an upwind approach was deemed inappropriate. Therefore, in the inter- 
est of simplicity, the predictor-corrector MacCormack scheme 19 was utilized for the present effort. 

MacCormack’s scheme can be applied to a scalar (or system of) conservation law(s) 


and is written as: 


where subscript ”i” and superscript 


du 

+ f£ =o 

(38) 

dt 

dX 


up 1 

= u" - AtVft 

(39a) 

■i(«r 

+ up 1 ) - \ AtAfl +l 

(39b) 


”n” denote spatial and temporal indices, respectively. Also, V 


and A denote fnst-oider backward and forward difference operators, respectively. In an analogous 


manner, we can rewrite the present system of equations as 


dq 

dt 


+ b' 


= 0 


(40a) 


where 


17 


b' = - L~ l b 


(40b) 


MacCormack’s scheme is then implemented as: 

q*+T = q'} + Atdqf ( 4 * a ) 

q- +l = \{<l n i + <7?^ + 2 Atd ^ (42) 

where, for example, 

dtfi = *>'? = - (L-'ty (43) 

In the predictor step, the vector function dq is evaluated using dependent variables computed at time 
level n and inverting the matrix L at each mesh point. This inversion is carried out using an efficient 
LU factorization 12 . Spatial derivatives in the vector b are approximated (conservatively) using 
first-onier backward differences where variables are again evaluated at the nth time level. Similar 
computations are performed during the corrector step, except that predicted values at time level 
n + i are used to perform the matrix inversion, and first-order forward spatial differences are used 
to approximate spatial derivatives. Because this is a central spatial difference scheme, additional 
numerical dissipation must be added to suppress unwanted oscillations. A simple fourth-order mod- 
el used by Warming and Beam 20 was used and implemented by modifying the corrector step Eq. 

(32b) above to give: 


where 20 


where CFL < 


«? +1 

—•|<N 

II 

+ qf*) + \At dqf^ + 

CsC<»d 4 q" 

(44a) 

a 4 ,; 

~ di + 2 

- 4 q n i + x + 6q? ~ *i?-i 

+ <fi- 2 

(44b) 



1 


(44c) 



C s g 





c,„ = 1 - CFL 2 


(44d) 


1 for stability. All solutions computed with this scheme were obtained with a CFL 


number of 0.9. 


d. Boundary Conditions 

For supersonic inflow and outflow, all dependent variables were specified and extrapolated, respec- 
tively. Conditions at subsonic inflow and outflow boundaries were treated by considenng the mvis- 
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cid and viscous equations separately. That is, for subsonic outflow, pressure was specified and densi- 
ty, velocity, and boundary-layer parameters were extrapolated. Mach number was then determined 
from velocity and the computed sonic speed. For subsonic inflow, the method proposed by Cooper 
et al 21 was used for the inviscid equations. By specifying inflow stagnation conditions, this method 
iteratively solves for the variables T e ,p e , and Ug using the equations (non-dimensional) 


y — 1 2 

T 0 ,e = Te + 


(45a) 


Pe = P Q,e 


C = U e - 


fef 


Pe 

Qe a e 


(45b) 

(45c) 


where Tis the temperature and a is the speed of sound. C used in Eq. (35c) is a “characteristic-like 
variable and is computed from information at the first mesh point inside the boundary 21 . Boundary- 
layer parameters 0 and T7 are specified and held fixed at the inflow boundary. 


e. Results 

The objective of this section is to present comparisons between computations obtained using the 
present interaction technique and measurements, as well as other computations. Subsonic and su- 
personic results are reported in separate sections where subsonic comparisons are all for steady dif- 
fuser flows, whereas supersonic computations are for both steady and unsteady channel flows. In 
all computations shown here, 51 equally spaced points were used in the axial direction. 


1. Subsonic Diffuser 

Axisymmetric subsonic diffuser flowfields investigated by Little et al 22 are compared with those 
computed using the present scheme (designated BL1D) in Figs. 5 and 6. The physical configuration 
consisted of several inlet pipe lengths (to give constant inlet boundary-layer thicknesses) and diffus- 
er half-angles, although only comparisons for the 12 degree, 21 inch configuration for inlet bound- 

* 

ary layer heights of 6 /R^ = 0.0034 (thinner inlet boundary layer) and <5 /R mlet = 0.0190 
(thicker inlet boundary layer) are reported here. It is almost embarrassing to report that computa- 
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tions resulting from the present interaction method required up to 30,000 time iterations to achieve 
convergence. However, it was suspected beforehand that this would be the case using an explicit 
numerical scheme for subsonic internal flows. Because the interaction methodology was of primary 
interest for the present effort (in particular, obtaining converged solutions which included bound- 
ary-layer separation), large iteration counts were considered an acceptable compromise between 
simplicity and efficiency. No attempt was made to optimize the code which was executed on a Sili- 
con Graphics, Inc. Personal IRIS 4d/30TG at a rate of 0.0016 cpu-seconds per time step per grid 
point. Therefore, a test case with 51 grid points requiring 10,000 time iterations resulted in approxi- 
mately 14 minutes of execution time. 

Also shown for comparison are computed parameters using a more classical interaction technique 
(herein designated as DUCFLO) where the inverse form of the steady integral boundary-layer equa- 
tions are iterated with edge velocity obtained from the constant mass-flow constraint. The inverse 
boundary-layer method used to obtain these results is that reported by Whitfield, et al 23 , although 
the DUCFLO interaction code (written by Whitfield), and findings generated by this code, have not 
been reported elsewhere. It should be noted further that this code can achieve converged solutions 
much more quickly than that using BL1D. However, the DUCFLO formulation is valid only for 
subsonic, steady flow and for this class of problems is generally the preferred technique with regard 
to computational resource requirements. 

Thinner Inlet Boundary Layer 

Figs. 5a-5e present comparisons between measured and computed distributions of static pressure, 
displacement thickness, momentum thickness, shape factor, and skin friction through the diffuser 
for the thinner inlet boundary layer case (integral lengths were formed using only kinematic proper- 
ties). Fig. 5a compares measured and computed static pressures (normalized by the inlet stagnation 
value), where the exit pressure (in BL1D) was adjusted until that at the inlet station matched the mea- 
sured value. Except for the region where diffuser diverence begins, computed pressures (from both 
BL1D and DUCFLO) are seen to compare favorably with those measured. It should be noted that 
the computed inlet pressure using BL1D was somewhat sensitive to the specified exit pressure. That 
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is, the exit pressure used to obtain the distribution shown in Fig. 5a was approximately 0.905, giving 
in inlet pressure of 0.60. Increasing the exit pressure to approximately 0.92 resulted in an inlet pres- 
sure of approximately 0.66. Therefore, cases reported in this section (using BL1D) were obtained 
by adjusting the exit pressure to match that at the inlet. 

Figs. 5b and 5c compare measured and computed “incompressible” displacement and momentum 
thickness distributions through the diffuser. It can be seen that both computational techniques over- 
predict and underpredict <Tand 0, respectively, although this agreement is considered reasonable. 
As a result of this over- and underprediction, computed shape factors are correspondingly high, as 
shown in Fig. 5d. It is reported in Ref. 22 that boundary-layer separation was not present in the 
experiment and none is predicted by the computations. This is illustrated in Fig. 5e which presents 
computed skin friction distributions (no measurements were available). However, exit shape factors 
in the range shown in Fig. 5d are an indication that considerable retardation in the velocity profile 
is present (i.e., the boundary layer is close to separation). This is illustrated in Fig. 5e which 
compares measured and computed (from BL1D) velocity profiles at the diffuser exit, where mea- 
sured profiles were obtained at three circumferential positions 120° apart. Although agreement be- 
tween measured and computed velocities is not particularly good at this axial location, considerable 
scatter exists in the data. Nonetheless, the computed profile is too thin and is also more retarded near 
the wall. However, as stated above, both measured and computed profiles are seen to be close to 

separation. 

Thicker Inlet Boundary Layer 

Fig. 6 presents comparisons between distributions of measured and computed (BL1D and DUC- 
FLO) diffuser parameters for the thicker inlet boundary-layer case. As stated previously, the exit 
pressure was adjusted until the computed inlet pressure approximated the measured value. Similar 
to the thinner inlet boundary-layer case, good agreement between measured and computed pressure 
distributions is indicated in Fig. 6a, although there is a “kink” in that computed by DUCFLO. This 
is apparently due to boundary-layer separation which is predicted by both BL1D and DUCFLO. 
It can be seen in Fig. 6b that good agreement between measured and computed displacement thick- 
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ness is obtained (except near the diffuser exit), although computed values are again slightly high. 
Similar comments can be made regarding measured and computed distributions of momentum 
thickness shown in Fig. 6c, although DUCFLO again exhibits a marked “kink” within the separated 
region which is not evident in the BL1D computation. As expected, computed shape factors are 
again too high corresponding to overpredicted displacement thickness. The extent of boundary-lay- 
er separation is shown in Fig. 6e which presents computed skin friction distributions for the thicker 
inlet boundary-layer case. The DUCFLO computations indicate a larger separated region than 
BL1D in that separation and reattachment occurs farther upstream and downstream, respectively, 
than does BL1D (again, no measurements were available). Figure 6f gives comparisons between 
measured and computed velocity profiles at the diffuser exit. Again, measurements were obtained 
at three circumferential positions around the diffuser exit. The computed velocity profile shows no 
reverse flow at this axial location (the computed boundary layer has reattached), whereas at least 
one set of these measurements indicate that the flow is separated. However, the agreement between 
measured and computed velocity profiles is considered reasonable. 

2. Supersonic Channel 

Steady and unsteady computations from BL1D for supersonic channel flows are compared to Navi- 
er-Stokes calculations in this section. Although supersonic nozzle flow calculations were reported 
previously 7 , the test case reported in Ref. 7 was for steady flow and compared only measured and 
computed wall static pressures. Attempts are made here to extend such comparisons to include the 
boundary layer, particularly for unsteady flow. 

Justification for comparing results ensuing from one computational technique to those of another 
comes from Hirsh 2 in reference to comments made in the Introduction; that is, computations result- 
ing from a Navier-Stokes analysis represent a higher level of approximation than those associated 
with the present methodology. The supposition here is that the technique employing the higher level 
of approximation is a better representation of the physics. While comparisons such as these are com- 
mon within the technical community, favorable agreement does not necessarily mean that results 
from the more approximate method represent reality; it just means that the two computations agree 
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with each other. The pitfall here is that while techniques utilizing a higher degree of approximation 
(i.e. more complete mathematics), may indeed represent more complete physics, the numerical 
method used to solve the resulting equations may be such that these physics are masked or otherwise 
lost Therefore, in keeping with previous comments regarding distinctions between physics, mathe- 
matics, and numerics, there can be no doubt that the Navier-Stokes equations are a more complete 
mathematical representation of the physics associated with supersonic internal channel flows. How- 
ever, we are assuming here that the numerical scheme used to give approximate solutions is yielding 
these physics to an acceptable level of accuracy which, of course, results in a better representation 

of the physics. 

The Navier-Stokes method used to obtain viscous solutions presented herein is that developed by 
Whitfield 24 and co-workers 25 " 26 . The particular version used most closely resembles that reported 
in Ref. 26 which has been modified to include explicit evaluation of viscous terms and extension 
of the solution algorithm to the so-called “modified two-pass scheme”. This code has as its basis 
an Euler solver which is an implicit finite-volume, formulation applying Roe s 27 approximate Rie- 
mann solver, and the higher-order extensions of Osher and Chakravarthy 28 , to compute the inviscid 
flux terms. The implicit operator is formed using Steger’s 29 flux vector splitting with the resulting 
system of equations inverted by application of Whitfield’s 30 two-pass or modified two-pass algo- 
rithm. The modified two-pass algorithm was applied in these computations. A brief description 
of the numerical scheme is included in the Appendix and the reader is encouraged to seek out the 
noted references for more details about the algorithm and the implementation. 

Comparison between BL1D and Navier-Stokes computations were made for two test configura- 
tions. The geometry analyzed was a transonic nozzle used as an AGARD 31 test case originally de- 
signed for evaluation of Navier-Stokes simulation capability relative to shock/boundary-layer in- 
teraction. To maintain isentropic, supersonic core flow, the nozzle pressure ratio was maintained 
below the second critical design pressure for the present computations. 
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Comparisons were made on this geometry for two test cases. The first comparison is made for steady 
supersonic flow subject to fixed inlet total pressure and total temperature and prescribed exit static 
pressure. The second case is for an unsteady flow and is created by linearly increasing the inlet total 
pressure as a function of time, whereas inlet total temperature is held constant and exit static pressure 
is again maintained near second critical design pressure. The geometry is shown in Fig. 7a. 

It should be pointed out that a problem arises when making comparisons between the Q ID analysis 
and the two-dimensional Navier-Stokes analysis. In the Q1D computation all flow-field parame- 
ters at a particular instant vary only as a function of axial location. However, each Navier-Stokes 
simulation produces a two-riimensional flow-field and determination of equivalent one^iimen- 
sional flow parameters for comparison is at best ambiguous. This results from the observation that 
determination of the boundary layer edge is not unique for a complex velocity profile. 


Steady Case 

For Navier-Stokes analysis, the AGARD transonic nozzle 31 was modeled using 153 X 30 mesh 
points in the axial and vertical directions, respectively. The grid spacing at the viscous wall, Ay/h 
is approximately .0002 ( where h is the channel half-height). This corresponds to a y + value of 
approximately 4. The Reynolds number based upon reference conditions and channel half-height 
is approximately 1.0 X 10 6 . The grid used and the boundary conditions specified in this simulation 

is shown in Fig. 7b. 

Identical boundary conditions and nozzle area variation were used to preform a corresponding simu- 
lation with the BL1D code beginning at axial location x=90mm, where inlet boundary conditions 
were taken from the Navier-Stokes simulation. Comparisons of computed distributions of Mach 
number, density, momentum thickness, and shape factor are shown in Figs. 8a-8d. Generally, agree- 
ment between computed core parameters from the two methods is considered good. However, mo 
mentum thickness and shape factor do not agree as well, where the largest disagreement occurs for 
150mm < x < 250mm. Fig. 8e shows velocity profiles in the nozzle region which illustrates that 
a unique definition of the boundary-layer edge is not possible thus causing the wide variations in 
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boundary-layer parameters computed from the Navier-S tokes results. This illustrates how the defi- 
nition of quantities such as momentum thickness and shape factor lose their significance in the con- 
text of complex velocity profiles such as those shown in Fig. 8e. 

Unsteady Case 

The AGARD nozzle geometry was modified for use as an unsteady test case. To produce a larger 
region within which comparisons could be made, the nozzle geometry was arbitrarily extended (in 
the axial direction) 4 nozzle heights as shown in Fig. 9. This geometry was then modeled with a 
185 X 30 grid similar to the grid used for the steady test case (Fig. 7b). Boundary conditions were 
as shown in Fig. 9. As stated previously, unsteady flow through the nozzle was initiated through 
temporal variation of the reservoir total pressure, shown in Fig. 10a. A uniform time step was se- 
lected such that the CFL number in the centerline region was near unity. This produced a maximum 
CFL number in the viscous layer in excess of 1 0 3 . This is expected to adversely affect temporal accu- 
racy of the simulation, but was deemed necessary in order to obtain results in a reasonable amount 
of CPU time. The impact of high CFLs occurring within large regions of the viscous flow field was 

not analyzed. 

A portion of the nozzle geometry (shown in Fig. 9) was analyzed with BL1D for comparison pur- 
poses. The Navier-Stokes simulation was used to define temporal boundary conditions at the inlet 
of the ”BL1D nozzle” (Fig. 9); these variations at x = 300 mm are shown in Fig. 10b. Comparisons 
between computed distributions of Mach number, density, momentum thickness, and shape factor 
atnon-dimensional times of approximately 500, 1000, 2000, and 3000 are shown in Figs. 11 through 
14, respectively. At all time levels, agreement between computed inviscid core parameters (a and 
b parts of Figs. 11 through 14) is considered very good. However, computed boundary-layer inte- 
gral lengths do not agree nearly as well, although the overall qualitative trends of boundary-layer 
behavior computed by BL1D are in good agreement with the Navier-Stokes results. 

The reader should note the “waves” and “wiggles” in the Navier-Stokes results shown in these fig- 
ures. As mentioned above, transforming two-dimensional results ensuing from the Navier-Stokes 
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analysis has proved to be somewhat challenging. These difficulties can be traced fundamentally to 
one’s definition of the boundary-layer “edge”; i.e., where the viscous region “ends” and the inviscid 
core “begins”. Integral lengths computed from the Navier-Stokes solution presented herein were 
generated by starting a search at the wall for the first maximum value of velocity (at a particular axial 
location and instant in time) which was then defined as the boundary-layer edge. Values of velocity, 
density, and Mach number at this y-location were then used to generate the various integrals. How- 
ever, significantly different values of edge quantities are obtained using another definition. For ex- 
ample, Fig. 15 compares computed momentum thickness distributions from the Navier-Stokes re- 
sults at the 1000 time-level (this corresponds to Fig. 12d) using two different edge condition 
definitions. For the second definition, the search discussed above was again performed to locate the 
first Umax in a particular profile. This value of velocity was then multiplied by 0.99 and another 
search conducted, again starting from the wall. The first index where the velocity exceeded the 
0.99iw value was then defined as the edge. This results in considerably different values of various 
edge quantities, the result of which is illustrated in Fig. 15. 

Therefore, although differences exist between computed unsteady core and boundary-layer quanti- 
ties resulting from Navier-Stokes and BL1D analyses, it appears that the Q1D approach is valid for 
this type configuration, at least for the conditions investigated. Again, however, additional efforts 
are warranted to investigate more consistent methods to interpret two-dimensional physics from a 
one-dimensional perspective. 
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III. PART 2 TRANSONIC FLOWS 

All previous discussion has been concerned with taking the system of equations (Eqs. ( 1) - (5)) and 
algebraically manipulating various terms in order to recast the original system in a form amenable 
for solution using relatively simple explicit numerical schemes. Although this approach has been 
shown to be capable of quickly yielding solutions of reasonable accuracy (at least for the cases pres- 
ented herein), the scheme possesses several distinct disadvantages. For example, many attempts to 
compute solutions containing near-discontinuous behavior of the dependent variables (i.e., shocks) 
with the explicit method have been unsuccessful, regardless of the type artificial dissipation model 
used (including that used in the following implicit scheme-see Section III. c). Because the ability 
to capture flows of this type is vital to any flow model which must operate in the transonic regime, 
it was evident that another approach must be pursued. 

Although not reported here, other explicit schemes have also been implemented, but again were not 
capable of capturing solutions with steep gradients. Therefore, the decision was made to implement 
an implicit scheme because of the inherent gains in stability bounds over those typically associated 
with explicit methods. It should be noted that this decision was originally prompted by the issue 
of numerical stability. It is recognized that if difficulties exist in obtaining numerical solutions for 
an ill-posed initial/boundary-value problem, it is probable that one’s choice of numerical scheme 
is not relevant. As stated previously, the present system of equations exhibit complex eigenvalues 
in flow regimes where shape factors (77) are high enough to induce boundary-layer separation. 
However, the magnitudes of the imaginary part of these complex eigenvalues are observed to be very 
small (one or two orders of magnitude less than the real pan). One could interpret this as meaning 
that the eigenvalue is “almost real” thus making the system “almost well-posed”. Nonetheless, as 
stated above, obtaining solutions using any scheme (especially those which are explicit) remains dif- 
ficult and therefore we cannot discard the possibility that the system is fundamentally ill-posed in 
certain flow regimes. Obviously, additional study is needed in this area. 

The algebraic complexities of the present system of equations limit the number of implicit schemes 
which can be used. For example, implicit schemes of Briley and McDonald 32 or Beam and Warm- 
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ing 33 can not be easily applied to this system because it cannot be written in fully conservative form. 
However, using a Newton-based implicit solver circumvents this limitation by discretizing the sys- 
tem in both time and space and then iterates the resulting system to convergence. The following 
discussion describes the Newton formulation as well as other issues resulting from its use. 

a. Formulation of the Newton Scheme 

Following the analysis presented by Whitfield 30 , a classical implementation of Newton’s method 
for finding roots to a nonlinear scalar function /(*) = 0 can be written as 

/ 'CO (x m+l - /*) = - fix" 1 ) ( 46a ) 

where m is an iteration parameter and 

f = # (46b) 

7 dx 

Now consider a system of nonlinear equations (each a function of several variables) written in very 
general form as 29 

F fx^Xfr •••■>x n ) = 0 
F ^ix j, ^2? ••• » Xfi) 0 

. (47) 

FniX\yX2i ^ 

If we consider F to be a vector function comprised of Fj , F 2 , .... F„, and * to be a vector function 
comprised of x } , x 2 , then the above system of equations can be written simply as 

F{x) = 0 

Newton’s method for such a vector F(x) can be written analogous to that for a scalar equation and 
solved as 

x m+ 1 = x m - [F '(Jt"*)]- 1 Fix”) (49) 

In the above equation, F \x) is the Jacobian matrix of the vector F(x) given by 
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where the elements of the Jacobian matrix are given by 

a ij (x) = -te- 


rn 


(51) 


That is, the (ij)* element of the Jacobian is given by the change in the i th element of the vector func- 
tion F(x) for a given change in the j* dependent variable. Because it is usually impractical to obtain 
the matrix inverse as the iteration proceeds, Newton’s method is usually implemented as 3 


F '(x m )(x m+l - x m ) = - F(x m ) ( 52 ) 

Now consider the vector function F(x) to instead be F(q), where the function F(q) is given by Eqs. 

( 1 ) _ (5). For example, the first element of F(q) is given by 


F, = A W) - Uef t (e<0 e ) + jdi^’ Rke) + - e ‘ u H - 0 (53) 

In implementing the implicit scheme, the formulation involving pressure (p e ) as a dependent van- 

able was used. Thus, the dependent variable vector is given by 


,2 C J- 


q = {Qe Pe ° ^ (54) 

The remaining elements of F(q) are defined by referring to Eqs. (2) - (5). The approach is to now 

discretize these functions using first-order backward temporal differences (implicit Euler) and se- 
cond-order central spatial differences, where the spatial differences are written at the (implicit) n+1 
time level. This results in 
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Of course, we assume that variables at time level n are known and we seek to solve for those at time 
level n+l. Therefore, within the framework of a Newton iteration, variables at the n time level are 
constant. However, each function F(q) above depends upon values of the dependent variable vector 
at spatial grid points i, i+1, i-1, all at the n+l time level. Thus, we can write 
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( 60 ) 


F(q n+1 ) = F{q^l q? +l , q^f) 

A Taylor’s series expansion of F (q n+1 ) in three variables results in the expression 



where 

A qn + \jn — qn + \/n+ 1 _ gn+l,m (63) 

Of course, the objective here is to perform sufficient iterations within a given time step such that 


which gives 


Aq n+lsn « o 


(64a) 


qii+ \/n+ 1 _ qn + ijn — qn+ 1 (64b) 

Eq. (62) is the Newton scheme used herein. It should perhaps be referred to as a discretized-Newton 
scheme because it is the discretized form of the equations to which the method is applied. The equa- 
tion is written at each interior mesh point which results in a system of block tridiagonal equations, 
where each block is a 5x5 matrix. A block tridiagonal solver written by Whitfield 16 is used to solve 
the system of equations. 

One should note that forming central differences on a non-uniform grid without the benefit of a cur- 
vilinear coordinate transformation has altered the formal second-order accuracy of the spatial 
discretization. This can be illustrated by considering a continuous function /(*) defined at discrete 
points Xi and forming Taylor expansions for /(x [+ / ) and fixi-i ) about /(.r,), or 
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For uniform spacing, the second term on the right-hand-side of the above vanishes and the usual 
central difference expression results. However, as evidenced by Eqs. (55) - (59), the present scheme 
simply uses the first term on the right-hand-side of the above and consequently incurs additional 
discretization error resulting from the non-uniform grid spacing. Fortunately, this additional “non- 
uniform spacing error” is typically smallest in regions where grid spacing is smallest, and larger else- 
where (at least using the stretching function discussed in Section in.d). Thererfore, the overall addi- 
tional error will be counter— acted by the second derivative term which will be small if grid points 
having the largest spacing have been placed in regions where the function is not changing rapidly; 
i.e., df/dx — 0. 

As mentioned above, only the interior points are updated using the Newton iteration. Boundary 
points are updated in an explicit manner at the conclusion of each complete time step in a manner 
appropriate with conditions at the boundary (i.e., subsonic inflow, supersonic outflow, etc.), as dis- 
cussed in Section n.d. 

The overall block structure is due to the appearance of the Jacobian matrices. Because of the algebra- 
ic complexity of the vector function F(q), the Jacobians are evaluated numerically. This procedure 
is discussed in the following section. 
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b. Computation of the Jacobian Matrices 

Consider the following evaluation of a Jacobian using the dependent variable vector at the i th mesh 
point: 

/V\ _ (65) 

W/i d (<h flifhflAfld] 

where 
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Therefore, we can write the (rs)* element of the Jacobian at the i* grid point (at the 
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m* iterate) as 
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The above relation states that at the i* axial location, the (rs)* element is computed by evaluating 
the change in the r* component of F(q) due to a given change in the s* component of q, holding 
q si+1 , q s i _j at their current m^-iterate values. 


The value of e used to compute the Jacobians as described above was approximately one-half of the 
reliable digits associated with the machine on which the code is executed, as suggested by Whit- 
field 34 . All solutions using the discretized-Newton scheme presented in this section were obtained 
on a Silicon Graphics, Inc. 4d/460 Power IRIS using double-precision floating-point operations. 

Therefore, for the present effort, £ ~ 10 " 6 . 


As might be anticipated, computation of the Jacobian matrices is rather expensive from an operation 
count point-of-view, particularly when considering that they appear within the Newton iteration; 
i.e., they should be recomputed at each m-iterate. However, it has been observed that convergence 
is not degraded significantly if the Jacobian updates occur infrequently. In fact, solutions presented 
here were obtained by updating the Jacobians about every five time-steps, where a time-step may 
include several m-iterations. This method of Jacobian updating drastically reduces the overall com- 
putational resource requirements. Note that this approximation only affects the convergence (i.e., 
to make Aq n+X/n ~ 0) at a particular time step and has no effect upon the time accuracy of the 
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scheme. This is because the right-hand-side of Eq. (62) is only affected by the discretization error 
associated with the temporal and spatial difference operators used in the original discretization pro- 
cess (in this case, first— order in time and second— order in space). 

c. Dissipation Model 

As stated above, second-order central differences were used for spatial discretization. As such, in- 
sufficient numerical damping is present which results in significant under- and overshoots in re- 
gions where the dependent variables exhibit near-discontinuous behavior (e.g., shocks). Therefore, 
additional artificial dissipation was needed to suppress these unwanted oscillations. 

Although several dissipation models have been tried (e.g., that used by Warming and Beam 20 ), the 
model proposed by Davis 35 has been the only one which yields the desired results. This is a “TVD- 
based” (total variation diminishing) model which determines the level of added dissipation using 
local flowfield conditions. The development of this model was motivated by the need to improve 
shock capturing capabilities of explicit schemes while preserving the simplicity of these methods. 
Much improved solutions were reported by Davis 35 in solving the two-dimensional Euler equations, 
and also by Causon 36 who used the model in three-dimensional inviscid flows, where the MacCor- 
mack scheme was used in both studies 35 - 36 . 

For the present discretized-Newton scheme, the additional dissipation was added to the dependent 
variable vector at the end of each Newton iteration (i.e., after each m-iterate). Although adding the 
dissipation at various other stages was tried, it was determined that adding it after each Newton itera- 
tion was the most robust and gave the overall best results. The model as implemented herein is given 
by the following step-by-step procedure: 

[q* + ^ jn + ^\ smoothed = [^ , " +1 ’ m ] unsmoothed ^ ^diss ^8) 

where 
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0(r) = min[mod(2r, 1)] = max(0, min(2r, 1)) 
where { a , b ) is a scalar product of a and b, and X is the maximum local eigenvalue. 

d. Results 

This section presents comparisons between measured 31 and computed results for steady, transonic 
converging-diverging planar channel flow for both symmetric and asymmetric configurations. In 
aHHirinn « results obtained by the present discretized-Newton scheme (hereafter designated 
BL1D-1), calculations ensuing from both Navier-Stokes (the same code previously discussed) and 
a Q1D Euler solver 37 are compared with the measurements. Solutions obtained with the discre- 
tized-Newton scheme were obtained using a Silicon Graphics, Inc. 4d/460 Power IRIS. Again, few 
attempts have been made to optimize the code for fast execution which proceeds at a rate of 0.01 8 
cpu-sec per grid-point per time-step. Calculations from BL 1D-I were obtained from an impulsive 
start by initializing the dependent variables to specified values and then marching in time to a steady- 

state. 

With regard to the comparisons which follow, the reader should note that. 
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1. X = 0 refers to the beginning of the convergent section. 

2. All lengths have been non-dimensionalized with h t , where h t = 0.1 meter 
for Case 1.2 and h t = 0.096 meter for Case 1.3. 

3. Static pressures are plotted as a ratio of local static to reference stagnation 
pressure (where p 0>oo = 96 kPA). 

4. Velocities are non-dimensionalized by the quantity u e , max , where 
Wejnax = 418 meters/second for Case 1.2 and u e , m ax = 403 meters/second or 

Case 1.3. 

Case 1 2 (Symmetric) 

A schematic of the experimental apparatus (which is symmetric) is shown in Fig. 16a. In attempting 
to computationally simulate the flowfield within a given physical configuration, it is important that 
all pertinent information about the experiment be given. Unfortunately, a complete description of 
the test apparatus was not provided in Ref. 31. In particular, the length of the section between the 
reservoir and the beginning of the converging section, as well as the distance between the diverging 
section and second throat, were not specified. Therefore, certain lengths had to be estimated where 
both estimated and actual values are shown in Fig. 16a. Inlet stagnation pressure and temperature 

were specified as p 0 x = 96 kPA and T O>00 = 300 K, respectively. These conditions resulted in a refer- 
ence unit Reynolds number of 5.345 x 10 5 per meter. Also, the reference length was defined to be 

R ref = 50 mm. 

For the present computations, 5 1 axial points were used where the points were clustered in the region 
of the shock. The clustering was accomplished using a hyperbolic tangent stretching function given 

by 

, _ tanh[<3(£ - 1)1 

* tanh(<5) 

where 


Here 6 is a scale factor that controls spacing and V max is the number of axial grid points. Use of 
such a stretching function allows a significant reduction in the number of grid points while maintain- 
ing good resolution in regions of high gradients. 
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Comparisons between measured and computed wall static pressure distributions are shown in Fig. 
16b. Symbols represent measured values whereas the solid line shows the computed distribution 
from the B11D-I code. Also shown for comparison is the computed pressure distribution using the 
inviscid method of Ref. 37. Both measurements and calculations indicate the presence of a shock 
of moderate strength. Similar to the subsonic calculations discussed in Section n, the exit static pres- 
sure in the computations was adjusted in order to match the measurements at some other axial loca- 
tion. In this case, the exit pressure was set such that the computed shock location was very near that 
indicated by the experimental data. Although both computational techniques show a reluctance to 
expand as quickly as do the measurements upstream of the shock, the inviscid calculations show a 
much higher rate of compression through and downstream of the shock than do those computed by 
the BL1D-I code. This “relaxing” effect is due to the presence of the boundary layer which, of 
course, the purely inviscid method cannot simulate. 

Shown in Fig. 16c is the computed skin friction distribution for this case. Whereas no experimental- 
ly determined values are available for comparison, Navier-Stokes results are shown and are seen 
to qualitatively agree with those computed by BL1D-I, although the separation points and the re- 
gions over which separation occurs are considerably different. Note that the computed distributions 
of Cyfrom both computational techniques are not smooth. In BL1D-I, this is due to similar behavior 

of the calculated “incompressible” shape factor 77 , which is shown in Fig. 16d along with the mea- 
sured distribution of this parameter. Although the computed distribution of H appears fairly smooth 
in Fig. 16d, a scaled-up plot shows significant “wiggles” exist upstream of the shock (many attempts 
to eliminate such oscillations have thus far been unsuccessful). The regions of experimentally and 
computationally determined reverse flow are shown in Fig. 16das the shaded and open areas, respec- 
tively. Further illustrations of this reverse flow are given in Fig. 16e which show measured and com- 
puted boundary-layer velocity profiles upstream and within the separated region, as well as at reat- 
tachment. Also shown for comparison are Navier-Stokes calculations for the same configuration. 
While the Navier-Stokes results show a clear superiority over those computed by BL 1D-I, the latter 
technique is seen to capture the overall trends inferred by the measurements. 
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Finally, Fig. 16f presents a comparison of measured and computed displacement surfaces within the 
perspective of the experimental test section. It should be noted that the vertical and horizontal scales 
of this figure are different which gives a distorted view of the minimum area region. It can be seen 
that computed results indicate a displacement surface considerably larger than that measured, al- 
though the overall qualitative trends associated with rapid boundary-layer growth downstream of 
the shock are reasonable. One should note that the computed displacement thickness within the sep- 
arated region grows to approximately 40% of the local channel half-height which implies that the 
viscous region occupies most of the channel (i.e., that the flow is fully developed). However, as 
shown by the velocity profiles in Fig. 16e, the rapid increase in displacement thickness can be attrib- 
uted to the fact that the velocity ratio ulu ^ is highly retarded over a significant distance from the wall, 
thus increasing the quantity 1-uJue used in computing the displacement thickness. Therefore, the 
flow should not be considered fully developed and that an inviscid core does exist (which can be 
verified from the velocity profiles shown in Fig. 16e). 

Case 1.3 (Asymmetric) 

The physical configuration for this case is similar to that given above except that the channel is now 
non-symmetric. Again, a complete description of the test apparatus was not provided, thus making 
it necessary to estimate the distance from the upstream reservoir to test section entrance. The exper- 
imental apparatus is presented in Fig. 17a which shows an airfoil-like blockage, but is on the lower 
wall only. Also shown in this figure is the computational domain used for the present calculations. 
In reality, this is a two-dimensional configuration which, of course, cannot be completely simulated 
using the Q1D approximation. However, because a large portion of internal flow configurations 
possess some asymmetries, it is of interest to compare results from non-symmetric measurements 
with those generated using the Q1D approximation in order to provide some measure of the validity 
with regard to using the Q1D approach, particularly for transonic flow. 

Shown in Fig. 17b are the actual and equivalent nozzle contours. The equivalent contour was gener- 
ated by using the actual area and deducing from this the required channel half-height distribution. 
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It can be seen that the effect of this “transformation” has been to lessen the wall slopes over the entire 
airfoil region, particularly at the aft end where the airfoil and tunnel wall coincide. 

Illustrated in Fig. 17c are comparisons between measured and computed static pressures where the 
experimental data shown were obtained on the lower wall (refertoFig. 17a) whereas those computed 
result from use of the equivalent channel as discussed above (Navier-Stokes results were not avail- 
able for this case). Also shown for comparison are wall pressures computed using the inviscid equa- 
tions 37 . Similar to the symmetric case, both computational schemes produce reasonable results up- 
stream of the shock, and again the inviscid computations indicate a more rapid compression 
downstream of the shock than do those associated with the present interaction scheme. Computed 
pressures from the present scheme again exhibit “wiggles”, but the oscillations seem to diminish 
axially downstream. It is important to note that the exit pressure for both inviscid and interaction 
calculations were adjusted such that the computed shock location approximated that observed on 
the lower wall. For both calculations, the use of an equivalent area distribution has the effect of un- 
der-expanding the flow in the region upstream of the shock. This is not surprising because the mea- 
sured flowfield “sees” a larger disturbance on the lower wall than that within the equivalent channel. 

Two-dimensional effects are shown in Fig. 17d which illustrates measured and computed pressures 
for both upper and lower walls. It can be seen that significant differences exist between measured 
upper and lower wall pressures and that the present calculations tend to represent an average of those 
measured, except with regard to shock location. This figure further amplifies the previous observa- 
tion that the equivalent channel (Q1D) approximation has the effect of causing a weaker shock (at 
approximately the same axial location as that observed on the lower wall). 

Similar to Case 1.2, a large region of separated flow was measured which is also predicted by the 
computations. Figure 17e shows computed skin friction which illustrates the region over which 
boundary-layer separation occurs in the calculations. The severity of the measured separation re- 
gion is further illustrated in Fig. 17f which presents measured and computed shape factor (W) dis- 
tributions. It can be seen that the calculated values of this parameter fall short of those measured 
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within the separated region. The reason for this underprediction is shown in Fig. 17g which illus- 
trates measured and computed velocity profiles as selected locations in the vicinity of the separated 
region. It can be seen that reverse flow regions associated with computed profiles are not as severe 
as those measured, thus accounting for the large shape factor. In addition, a correspondingly large 
displacement surface compared to the complications is shown in Fig. 17h. Whereas the computed 
displacement surface downstream of the shock is less than that measured, the calculations are con- 
sidered reasonable. This is particularly true when one considers the level of approximation used in 
the computational method. Similar to Case 1 .2, this case does not appear to be fully developed thus 
allowing the present interaction method to be applicable, at least for these flow conditions. 
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IV. SUMMARY AND CONCLUSIONS 

Viscous-inviscid interacting internal flowfields have been computed for subsonic, transonic, and 
supersonic quasi-one-dimensional (Q1D) configurations. Viscous-inviscid interaction was 
achieved by directly coupling the unsteady Q1D Euler equations with integral boundary-layer equa- 
tions for unsteady, two-dimensional, turbulent flow over impermeable, adiabatic walls. The cou- 
pling methodology differs from that used in the past in that the coupling is carried out directly 
through the equations as opposed to solutions of the different equation sets. Numerical solutions 
to the coupled system of equations were obtained using the explicit MacCormack scheme as well 
as an implicit discretized-Newton scheme. Computed solutions have been compared with measure- 
ments as well as Navier-Stokes and inverse boundary-layer methods. Although differences exist 
between measurements and solutions computed from the Q1D approach, and also between those 
computed by other methods, overall qualitative trends for both steady and unsteady test cases are 
predicted with reasonable reliability. 

An analysis of the eigenvalues of the coefficient matrix associated with the quasi-linear form of the 
coupled system of equations used in the Q1D approach indicates the presence of complex eigenva- 
lues for certain flow conditions (in particular, values of shape factor, TT, high enough to cause bound- 
ary-layer separation) thus allowing exponential growth in the solution(s). Although reasonable 
solutions are obtained numerically, it is believed these complex eigenvalues contribute to the overall 
difficulty observed in obtaining numerical solutions to the coupled system of equations. It is further 
postulated that these complex eigenvalues are a result of empirical and analytical approximations 
used in the integral boundary-layer technique. 

The original formulation of the system of equations 7 (for use in the explicit scheme) made the as- 
sumption that stagnation temperature was constant from the upstream Teservoir reference state to 
that associated with local conditions. However, this assumption was shown to be inconsistent within 
the context of unsteady flow, even for the case of no heat transfer. As a result, new formulations 
of the equations were derived which do not depend upon an assumption of constant stagnation tem- 
perature. It was shown that the degree of error associated with the original formulation was small. 
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at least for the case considered. However, this error will increase in direct proportion with the degree 
of unsteadiness associated with a particular flowfield. 

Comparisons between computations using the Q1D approach and those resulting from a Navier- 
Stokes analysis show generally good agreement with regard to the overall qualitative aspects of the 
flowfield. Some of the observed discrepancies can be attributed to improper interpretation of the 
two-dimensional Navier-Stokes results from the perspective of quasi-onedimensional physics. 

In particular, locating the boundary-layer edge in a Navier-Stokes generated velocity profile proved 
to be the “Achilles heal” of the comparison process. Nevertheless, results generated from the Q1D 
approach for both steady and unsteady supersonic channel flows are considered good enough to pro- 
vide preliminary estimates of internal flowfield behavior, at least for configurations of the type con- 
sidered herein. 

Implementation of an implicit discretiaed-Newton scheme for numerically solving the coupled sys- 
tem of equations was originally motivated from the standpoint of numerical stability. This was 
deemed necessary for the method to be capable of addressing transonic flow; i.e„ flows with shocks. 
However, as stated above, appearance of complex eigenvalues for the coupled system gives rise to 
possible exponential solution growth due to the equations being ill-posed for solution in time as an 
initial/boundary-value problem, at least in certain flow regimes. In spite of this, reasonable steady- 
state solutions have been obtained, even for highly separated flows (i.e., for very high values of H). 
It is concluded that even though the equations appear to be ill-posed in these regions, growth of error 
is slow enough to allow reasonable solutions. 

It is acknowledged that “there’s nothing more dangerous than answers that look about right,” and 
that more study is warranted with regard to both unsteady and multi-dimensional flow, simulation 
using the Q1D approach. However, it is felt that the Q1D approach is valid and the tools generated 
by this effort are useful in the regimes which have been investigated, particularly when considering 
the computational resources required to generate results from other methods. For example, the Nav- 
ier-Stokes code used in these studies executes at approximately 6.0 x 10~ 5 cpu-sec per time-step 
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per grid point on a CRAY Y-MP. Therefore, for 3,000 time steps and the grid used (185 x 30 x 2), 
this resulted in approximately 30 cpu minutes on this machine. As stated earlier, all QlD-based 
codes were executed on engineering workstation-based hardware and required much less computa- 
tional time to execute. As stated in the Introduction, whether or not the Q1D approach is appropriate 
is greatly problem dependent, as evidenced from the steady-state nozzle results which was inherent- 
ly two-dimensional within the nozzle contraction/expansion region. If preliminary, engineering ac- 
curacy results will suffice, results shown herein indicate that this can be achieved using this ap- 
proach. 
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Figure 3. Computed Unsteady Flow Parameters (All F ormulations) 
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Figure 3. (Continued) 

(c). Stagnation Temperature Variation 
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eigenvalues 


Figure 4. Eigenvalues ofL N (lines denote eigenvalues of 
complete system, symbols denote u e , u e +a e ; u e -a e ) 


(a). H = 1.2, 6 = 0.001 



(b). H = 1.5, 0 = 0.001 
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eigenvalues 


Figure 4. (Continued) 
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Figure 4. (Continued) 
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Figure 4. ( Continued ) 

(g). h = 1.5, e = o.oi 
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Figure 4. (Continued) 
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Figure 4 . (Continued) 
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Figure 5. Measured and Computed Subsonic Diffuser Parametei s 

' y (Thinner Inlet Boundary Layer) 


(a). Wall Static Pressure 
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Figure 5. (Continued) 

(c). Incompressible Momentum Thickness 
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(d). Incompressible Shape Factor 
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Figure 5. ( Continued) 

(e). Skin Friction 
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Figure 6. Measured and Computed Subsonic Diffuser Parameters 

6 (Thicker Inlet Boundary Layer) 


(a). Wall Sialic Pressure 
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Figure 6. (Continued) 


(c). Incompressible Momentum Thickness 
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Figure 6. (Continued) 

(e). Skin Friction 
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(f). Exit Velocity Profiles 
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Figure 7 . AGARD Geometry 

(a). Test Apparatus 
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Maclt Number, M 


Figure 8. Navier-Stokes andBLID Comparisons (Steady Case) 
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Figure 8. (Continued) 


(c). Momentum Thickness 



x (mm) 


(d). Shape Factor 



BL1D 


t' o Navier-Stokes 

1.8 r 


1.7 r 



x (mm) 


67 


(imu) X 


Figure 8. (Continued) 

(e). Velocity Profiles in Nozzle Region 
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Figure 9. Computational Grid and Boundary Conditions for 
Modified AGARD Geometry ( Navier-Stokes — Unsteady) 
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Various Non-Dimensional Parameters 


Figure 10. Time Variance of Navier -Stokes Parameters 
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Figure 11. Navier-Stokes and BL1D Channel Parameters at t - 500 

(a). Mach Number 
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Figure 11. (Continued) 

(c). Momentum Thickness 
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Figure 12. Navier-Stokes and BL1D Channel Parameters at t - 1000 
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Figure 12. (Continued) 


(c). Momentum Thickness 
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Figure 13. Navier-Stokes andBLID Channel Parameters at t = 2000 
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Figure 13. (Continued) 


(c). Momentum Thickness 
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Figure 14. Navier-Stokes and BL1D Channel Parameters att = 3000 

(a). Mach Number 


1.62 
1.61 
1.60 
1.59 
1.58 

^ A57 j- 

[ 

1.56 V 

I 

1.55 Y 

L 

1.54 | 

l 

1.53 i - 


1.52 


— computed (B LID) 
o computed (Navier-Stokes) 


'Oq 


° o o 



9-e- 


o o o o ' 


300 350 400 450 500 550 600 650 700 750 

x ( mm) 

(b). Density 


0.77 
0.76 
0.75 
0.74 
£ 0.73 
0.72 
0.71 
0.70 
0.69 


0.68 


computed (BL1D) 
computed (Navier-Stokes) 


300 


350 


400 


450 


500 550 

x (mm) 


600 


650 


700 


750 


11 



Figure 14. (Continued) 


(c). Momentum Thickness 



(d). Shape Factor 
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Figure 15. Momentum Thickness Distributions at tZ 1000 

(Navier-Stokes) 
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Figure 16.AGARD Test Case 1.2 

(a). Experimental Test Configuration 



(b). Measured and Computed Wall Pressure 
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Figure 16. (Continued) 

(c). Computed Skin Friction 



x/h ( 

(d). Measured and Computed Boundary-Layer Shape Factor 



* 


81 


(mm) y (mm) 


Figure 16. (Continued) 

(e). Measured and Computed Velocity Profiles at Various Axial 

Locations 
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Figure 16. (Continued) 

(f). Measured and Computed Displacement Surfaces 
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Figure 17. AGARD Test Case 1.3 

(a.) Experimental Test Configuration 



(b). Equivalent Wall Height Distribution 
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Figure 1 7. (Continued) 

:). Measured and Computed Wall Pressures (Measurements from Lower Wall Only) 



xlh t 

(d). Measured and Computed Wall Pressures (Measurements from Lower Wall Only) 




Figure 1 7. (Continued) 

(e). Computed Skin Friction 
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Figure 17. (Continued) 

(g). Measured and Computed Velocity Profiles at Various Axial Locations 
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Figure 17. (Continued) 

(It). Measured and Computed Displacement Surfaces 
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APPENDIX 


I. Elements of the L and N Matrices 

Recall that the system of equations can be written as 


or, in quasi-linear form, 


T dq _ . 

L Tt ~ b 




(A.I.l) 


(A. 1.2) 


Elements of the L and N matrices are given below for all formulations and were obtained using the 
auxiliary relations reported elsewhere in this Appendix. Each particular element is referred to by 
either lower case € (for L) or n (for N) using double subscripts (ij), where i and j represent the specific 
matrix row and column, respectively. 

Listed first are the L and N matrix elements resulting from the dependent variable vector given by 

q = ( g e u e M e 6 TJ) T 

Although the ratio T 0 e = f 0e /f 0> » appears in several terms, all computations using this formula- 
tion have been made using the assumption T 0 , e = L It should be further noted that additional terms 
involving derivatives of this quantity were neglected in the original formulation, i.e., 


dT 


0,e 


dT, 


d t 


0,e 

dx 


= 0 


I Matrix Using M. as P endent Variable (T<y assumed constant) 

l u = U<0 (H & . — Hg) 

1 1 2 = Q^^d’ 

/ 13 = iQelleMJ) ( C 5 H + C 6 ) 

^14 = Qe^eifiy ~ Hq^) 
l l5 = Q e uJ)(1 + C 5 Ml) 

1 21 = 6u 2 e ( 1 + ~ Hg e ) 

1 2 2 = 2Q e U^l + H d . - 

^23 = 2 Q e ulMed(C 5 H + C 6 ) 
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h\ - Qeutil + H d . - H Qq ) 

l 25 = Q e uld{ 1 + C 5 Mj) 

hi = A 

hi = 0 

/ 33 = - iQeMedw (C X H + C 2 ) 

1 34 = Qe w Hd' 

1 35 = - Q<8w( \ + C X M\) 

/ 4 1 = M 

hi ~ QeA 

/ 43 = - lQ e u e MJ)w ( c x H 4 - c 2 ) 

I44 = - QeU e wH d . 

/ 45 = - QeU^w(\ + C X M \ ) 

^51 = A^ieTft'g + 2 U «) 

^52 = Qe^eA 

l 53 = - 2{E e + Pe)v$M e ( c x H + c 2 ) 

^54 = (E e Pe)wHfi* 

l 55 = - (E e + Pe)wd (1 + c x M 2 ) 


N Matrix Using M» as Dependent Variable (Tn, 

n n = du 2 e 

n l2 = QeU^(H 6 . + 2) 

n 13 = 0 
H 14 = Qe^e 

n 15 = 0 
fi 2x = Oh^H q* 

n 22 = Q e uld{^>H Q . + 2 H y - 2yj 

3 

n 2 3 = QeUeQjM^ 

^24 = QeUeHfp 



AQfM 6 Tq e 


assumed constant ) 
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_ dHg* 
n 25 = QeUe6 ~W 

rt 31 = UeA 

IX 32 = Qe& 

n 33 = - 2QeUeM^W(c x TT 4- c 2 ) 

« 34 = - Q e u e wH 6 . 

n 35 = - QeUjam, 1 + c x Ml) 

n 41 = A[(y - D^o.* + 

71a2 = 

__ y — \Ag e M err 

n 43 = - 2Q e U 2 e dM e w{C x H + c 2 ) y y2 y o,e 
«44 = - QeU 2 e wH d . 

n 45 = - Pe«e^vv(l + 
n 51 = UeA^ye e T 0e + \ u i) 
n 52 = A(Ee + Pe + G^e) 

_ s QeUMArj, 

n 52 = - 2(£ e + p e )u e w6M e (CiH + c 2 ) y| 7 o,e 

« 54 = — (Ee + Pe)UeWHd' 

n 55 ~ “ ( £ e + />e)H e H#(l + CjAfl) 


Also, the following definitions have been used. 


Cl = 0. 113 
c 2 = 0.290 
c 3 = 0. 185 
c 4 = 0 . 150 

c 5 = c\ ~ C 3 


c 6 = c 2 - c 4 


A = w (R - 6*) 


w = 



2 , planar 

• 1 ( 7 ? — 6*), axisymmetric 
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w = 


{ 


w, planar 
2 w, axisymmetric 


~ = — 

e * ~ 6e °' e ~ y(y ~ We y(y ~ D M 2 e Y<y ~ D 

E e = Qe(e e T 0<e + = Qe(^-t ' ^ 


1) 


+ i u * 


_ QeTe 
p e v 


E e + p e = Qe{ye e + 
y - 1 


fe ~ 1 + 


Ml 


Listed next are th e L and N matrix elements using Mach number as a dependent variable, where no 
assumption regarding T 0 ^ = constant is made. For convenience, most elements are written in terms 
of those given previously. Matrix elements given below are written with superscript M to indicate 
that Mach number is used as the dependent variable and that stagnation temperature is not assumed 


constant. 


L Matrix Elements Using Mach Number as Dependent Variable 
(To >e not assumed constant) 


iM - I 

Ml " ‘11 

/M _ / 

12 ~ ‘12 

tM — } 
M 3 ‘13 




1 14 “ * 14 

iM _ / 
*15 “ *15 


l M = / 
*21 * 


21 


iM _ / 
*22 ” *22 

iM 

*23 


— / 

fo* — <23 


_ / 

‘24 “ *24 

jM - j 
*25 “ *25 

iM _ / 
*31 “ *31 
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/M . 
*32 ' 

= hi 

/M . 
*33 ' 

- hi 

*34 

= / 34 

/A/ 

*35 

»r> 

ro 

ll 

iM 

*41 

“ ^41 

,Af 

*42 

= / 42 

;Af 

*43 

= / 43 

M 

*44 

= I44 

M 

*45 

= ^45 

M 

*51 

= Ak h 

iM 

*52 

= 2A/ 

M 

*53 

= — j 

M 

*54 

II 

iM 

*55 

= hs 


2 AQeU % 


2 (E e + Pe)wQM e {c x H + c 2 ) 


N Matrix Elements Using Mach Number as Dependent Variable 
(To,e not assumed constant) 

_ „ 

«H _ "ll 
«12 = " 12 
" l3 — " l3 

M A/ — y% 
n \A ~ n 14 

«15 = "15 
n 2\ = "21 

"22 = "22 

«23 = "23 

„M _ „ 

«24 ~ "24 

„M _ „ 

"25 ~ ”25 
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*31 — ”31 
*32 = ”32 
*33 = ”33 
”34 = ”34 
”35 = ”35 


u 0 Au l 
yjM — All -f* ■■ 

”41 AUe + vM2 


= 2 g e UeA + 


yM'l 

'2 *Aq gUg 


yM\ 


2 # ^ 

”^ 3 = “ “ ZQeulwQMeicfl + C 2 ) 


‘43 

„W _ . 

*44 " ”44 

*45 = ”45 

n?l = 


”52 = A(£ e + Pe + 2QeUekM) 


(y - DM 3 


«53 = - 7^7tL - 2 (^ + Pe)u e yv6M e {c x H + C 2 ) 

*54 = ”54 

„M _ „ 

*55 “ ”55 


Note the following relations have been used: 


~ ~ + 


2 y(y - 1)A/ 3 

*'« - l + ^ijsf - i (1 - *> + *« 

Finally, the L and N matrix elements resulting from the dependent variable vector given by 

Q = {Q e Me pe @ 

are listed. As previously discussed, no assumptions of constant stagnation temperature are required. 
Matrix elements listed below are written with superscript p to signify that pressure was used as a 
dependent variable. Again for convenience, most elements are written in terms of those previously 

listed. 
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L Matrix Elements Using Pressure as Dependent Variable (T„, not assumed constant) 


[P = 
11 

l P = 
*12 

IP = 

13 

IP = 

14 

IP = 

15 

IP = 
21 

IP = 
22 

]P = 
23 


, + ii M* 

'll + 2 13 Qe 

, Me 
hi + hs~uf 

1 / 

“ 2 13 Pe 


14 


hs 

z +1, Ms 

*21 + 2 23 Qe 
, We 

^22 + '23 M e 


_I; 

2 23 Pe 


/P 

24 

/P 

25 

IP 

31 


= / 


24 

'25 

1 , W* 
'31 + 2 /33 e7 


/p 

32 


, M e 

*33177 


/P 

33 


1, ^ 
- 2*33 


II 

Tf 
cx m 

l 24 




11 

yn 

txn 

l 35 





h\ 


1 . 

M, 

ip = 

41 

+ 

2 43 Qe 


hi 


'43 

M e 

IP = 
42 

+ 

■uf 


IP 

43 


_ _I; MS 

2 43 Pe 


/P 

44 

IP 

45 


- ^44 

= '45 
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1 rut ** 

^5i = l Au ‘ ~ ~oi^ Ee + Pe)w6 ( c ' H + Cz) 

M 2 — 

IP = QelieA. - 2-^{E e + p e )w6(c { H + c 2 ) 
jZ * 


4 Ml 


IP 3 = + pf (£e + Pe)mc { H + c 2 ) 


*54 “ *54 
?5 = '55 


N Matrix Elements Using Pressure as Dependent Variable (To, e not assumed constant) 

"it = "n 
"ia = "i2 
"l3 = "13 
* P 14 = "14 
"iS = "15 

„ 1 M e 

n 2\ ~ "21 + 2 " 23 Qe 

n Me 

"22 = "22 + "231Z7 
. 1 M e 

"23 “ “ 2" 23 


"24 = "24 
"25 = "25 

„ 1 M e 

n 31 = "31 + 2"33?7 

„ M e 

«3 2 = "32 + "33 TT e 

p \ Me 

"33 ” “ 2"33 J 7 

nP M = "34 
"35 = "35 

n p 41 = Aul - ulMldncJJ + c 2 ) 
n p 42 = 2QSeA - 2Q<UeMl6w{c x TJ + c 2 ) 
n p A3 = A + ^jf-M 2 e ew{c x H + c 2 ) 
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nP U = n 44 
n 45 = * 45 

n 51 = + Pe)Wd(ClH + Cl) 

nP 52 = A{QeU 2 e + E e + p e ) ~ 2M 2 (E e + Pe)wQ{c-JJ + c 2 ) 

n P = -t-rrAUe + ^M 2 (E e + Pe^ic^ + C 2 ) 

"53 y - 1 P* 

n 54 = "54 

"55 = ” 55 

IL Auxiliary Relations 

H y = (1 + c x M 2 )H + c^M 2 
Hq 9 = ( c 3 /7 + c 4 )Mi 

1 _i[io 7 - 77 ' - 1 

Wa#.- 0 = 1 • 48061 + 3 • 83781e_2 ^ + 0 • 33 “ 8 . 5484 tan [ 1-23 

- (0. 33 - yf-fjtanh^l . 2874 x l<r 6 )(10 7 - ff >’ 45761 ] 


H e - = 


(H-Om.-o + 0 ■ 028M « 

1 + 0.014M 7 


0.3 


— 1.337T 


+ (1 . 1 x 10 4 ) 


- <log 10 J& 9 )'-’ 4+0 ' 3,n ' 


| tar| h ( 4 0 .^ 75 ) 


- 1 


8 - 1 - 0 • 92M i-tanh[l . 49(// -0.9)] 

0 7.09 + M 2 


2 = 
2 


- 0 . 01 167 e 


-0.03877 3 + 0.0115 + W + < 8 9 - 0X l0 ' 8yl 


.60377 


(1 + 0.025A/1- 4 ) 


Note that in the above relation, 


ACFD = mReJ 

m = 65077 - 743 
n = - 1 . 5977 + 


} ” 

1 .45 y 


< 1.6 
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ACFD = me 


m 

n 


nT^e^ 


3 . 25e 0 M5TJ * 


77 

10000 


0. 0017 


} 


77> 1 .6 


where 


= 5 = A 

Re, c f Je 


in. Details of Navier-Stokes Computational Procedure 

The Navier-Stokes simulations of the AG ARD transonic nozzle geometry were made by application 
of the Euler solver developed by Whitfield 24 and Arabshahi 26 modified to include viscous effects. 
Although this code was developed as a full Navier-Stokes code, stream wise viscous terms were ne 
glected in these two-dimensional simulations. 

The basic Euler solver is an implicit, finite-volume, formulation applying Roe’s 27 approximate Rie- 
mann solver, and the higher-order extensions of Osher and Chakravarthy 28 , to compute the inviscid 
flux terms. The implicit operator is formed using Steger's» flux vector splitting with the resulting 
system of equations inverted by application of Whitfield’s* two-pass or modified twev-pass algo- 
rithm. The modified two-pass algorithm was applied in this computation. 

Viscous fluxes were added as an explicit source term patterned after the implementation in the 
PARC3D 21 Navier-Stokes code. The viscous portion of the Reynolds-averaged Navier-Stokes 
equations can be expressed in nondimensional conservation law form as 

j \j2i 

Re dtj 

As usual, the Reynolds number. Re, is defined by reference sound speed, length, density, and vis- 
cousity. The viscous flux vectors are defined as 


G j 


= J 


0 

r ij 

« k r jk ~ 4j 
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where /is the Jacobian of the coordinate transformation. The viscous stress tensor and the heat flux 
vector are defined as 

d£j dUi dUk 1 ■■ dU m 

T V — ^ dx k d.Xf, dXj J dX; dx m 

K d JhdT_ 

q J (y - 1 )P r dx k dx k 

The viscous flux is then approximated at each cell face in a series of directional sweeps with the result 
summed into the residual computed for the original Euler algorithm. 

The turbulence model applied is also based upon the algebraic model implemented in the PARC3D 
Navier-Stokes code. In wall-bounded regions of the flow field, a Baldwin-Lomax 38 algebraic 
model of turbulent viscosity is applied. In regions not bounded by a solid surface, a vorticity based 
model, as developedby Thomas 39 , is applied. In this algebraic model, viscosity and thermal conduc- 
tivity are modified as a function of turbulent viscosity as 

P total = A* + 

K total K . f*T 

Pr Pr P rT 

ptj = R qq1 2 oj 

where PT is the turbulent viscosity, P rT is the turbulent Prandtl number, I is a turbulent length scale 
and to is the local vorticity. 
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NOMENCLATURE 


a speed of sound 

A area 

b right-hand-side vector 

Cf skin friction 

c v constant volume specific heat 

C function used in dissipation model; also used as constant for inflow boundary conditions 
D dissipation integral 

e internal energy 

E total energy 

f scalar function 

F vector function 

G function used in dissipation model; also denotes viscous flux vector 

h channel half-height 

H “incompressible” shape factor <5 jd 

Hy shape factor, <3 */6 

Hq 9 shape factor, 0 Q /0 

Hq. shape factor, 6* /0 

K thermal conductivity 

L temporal derivative coefficient matrix 

€ ~ (ij)* element of L matrix 

M Mach number 

N spatial derivative coefficient matrix 

N m ax total number of grid points 

n ij (ij)* element of N matrix 

p pressure 

Pr Prandtl number 

q dependent variable vector, also denotes heat flux vector in Navier-Stokes code 
r function used in dissipation model 

R radius; also used to denote perfect gas constant 
Re Reynolds number 

t time coordinate 

T temperature 

u velocity 

u T friction velocity ( u? = t w/q) 
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x axial coordinate; also used to denote dependent variable in Section IH.a 
y coordinate normal to wall 

( , oyuA 
y = —p~J 

y ratio of specific heats 

<5* displacement thickness 

d T “incompressible” displacement thickness 

< 5 ; V 

A forward difference operator 

V backward difference operator 

6 momentum thickness 

6 “incompressible” momentum thickness 

6 * energy thickness 

6 e density thickness 

\ eigenvalue; also denotes second coefficient of viscosity 
p absolute viscosity 

v CFL number 

| parameter used in grid stretching; also denotes curvilinear coordinate used in Navier-Stokes 
code 

q density 

t shear stress 

<J> function used in dissipation model 

to vorticity 

Subscripts 

e boundary-layer “edge”, or inviscid core value 

i axial index 

m iteration parameter 

max denotes a maximum value 

n denotes the n* equation or dependent variable; also denotes time level 
ref reference condition 

T denotes a turbulent quantity 

0 stagnation condition 

w wall 

oo reference condition 

Superscripts 

k 0 for planar flow, 1 for axisymmetric 
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M signifies Mach number used as dependent variable and To >e * constant 
p signifies pressure used as dependent variable and To >e * constant 
A dimensional quantity 
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